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We compute the virial coefficients, the contact parameters, and the momentum distribution of a 
strongly interacting three-dimensional Bose gas by means of a virial expansion up to third order in 
the fugacity, which takes into account three-body correlations exactly. Our results characterize the 
nondegenerate regime of the interacting Bose gas, where the thermal wavelength is smaller than the 
interparticle spacing but the scattering length may be arbitrarily large. We observe a rapid variation 
of the third virial coefficient as the scattering length is tuned across the three-atom and the atom- 
dimer thresholds. The momentum distribution at unitarity displays a universal high-momentum 
tail with a log-periodic momentum dependence, which is a direct signature of Efimov physics. We 
provide a quantitative description of the momentum distribution at high momentum as measured by 
P. Makotyn et al. [Nat. Phys. 10, 116 (2014)], and our calculations indicate that the lowest trimer 
state might not be occupied in the experiment. Our results allow for a spectroscopy of Efimov states 


in the unitary limit. 

PACS numbers: 67.85.-d, 67.10.-j, 67.10.Hk, 34.50.Cx 

I. INTRODUCTION 

The classical problem in Newtonian mechanics of find¬ 
ing analytical stable orbits of three planets has only a 
few known solutions which, moreover, exist only for cer¬ 
tain special conditions pQ. Curiously, a generic analog 
of this problem in quantum mechanics, three neutral 
bosons of mass m with a universal short-range inter¬ 
action, was solved analytically by Efimov four decades 
ago [2H1] • If the effective range of the interparticle inter¬ 
action (set by the van der Waals length £ v dw ) is small 
compared to other length scales, the two-body interac¬ 
tion is solely characterized by the scattering length a. 

In the unitary limit of infinite scattering length, Efi¬ 
mov found an infinite number of three-particle bound 

states with energy = — (e ,r / So ) _2j , where re* is 

a universal three-body parameter, j is an integer, and 
e 7r/so 22.7. Originally predicted in nuclear physics, 
Efimov states were observed in atom- loss experiments 
and radio-frequency spectroscopy measurements in Bose 
gases PI, and three-component Fermi gases EMU , as 
well as mass-imbalanced mixtures [Mm. Recently, the 
Efimov trimer of 4 He has also been observed experimen¬ 
tally [IS. Efimov physics is thus a very general phe¬ 
nomenon, and in addition to these many experimental 
realizations, Efimov states are also predicted to exist, for 
example, in p -wave interacting quantum gases tm, in 
condensed-matter quantum magnets [ 20 ] . and in mass- 
imbalanced two-component Fermi gases [2D H2] Most 
theoretical work focuses on few-body aspects of Efimov 
physics and experiments are commonly explained using 
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few-body theory [J] . In this paper, by contrast, we study 
the interacting Bose gas at finite density and temperature 
and establish signatures of three-body correlations in this 
system. The aim of this work is to provide a first exact 
quantitative calculation of three-body effects on a many- 
body system, which is relevant from a fundamental point 
of view and can be compared with current experimen¬ 
tal work. The main result of this paper is a calculation 
of the full momentum distribution in the nondegenerate 
limit for arbitrary values of scattering length and three- 
body parameter re*. 

Most experiments on Bose gases in equilibrium are re¬ 
stricted to the weak-interaction regime with few excep¬ 
tions g3H2U. This is due to enhanced three-body losses 
at finite density with increasing interaction strength 
which deplete the gas with a rate N = L 3 n 2 N [?] . At zero 
temperature, the loss coefficient scales as L^,{T = 0 ) ~ 
a 4 [25M27] and saturates to L 3 (T) ~ ~ 1/T 2 [2BJ in 

the unitary limit at high temperature (At = \J2n/mT 
denotes the thermal wavelength, and we set h = k B = 1 ). 
In a series of recent hallmark experiments, the loss rate 
of a strongly interacting Bose gas at finite temperature 
was measured [22H52| . It turns out that in the nonde¬ 
generate limit nX B <C 1 , the two-body scattering rate 
72 = nav ~ uXt is much larger compared to the three- 
body loss rate 73 = L 3 n 2 ~ n 2 X B , so that the gas can 
reach an equilibrium state before a significant fraction 
of particles is lost. Indeed, this is corroborated by a re¬ 
cent experiment that quenches a weakly interacting Bose- 
Einstein condensate to the unitary limit and measures 
the momentum distribution, which approaches a station¬ 
ary equilibrium distribution shortly after the quench m 
In the nondegenerate regime, we furthermore expect that 
the system is thermodynamically stable . The ex¬ 

periments [ 2 !M? 2 j demonstrate that strongly interacting 
Bose gases, for which three-body correlations turn out to 
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be essential, can be experimentally prepared and studied. 

Given this experimental progress, there is a need for 
an accurate quantitative description of a strongly inter¬ 
acting Bose gas which includes three-body effects. Yet, 
typical many-body methods, such as the T-matrix ap¬ 
proximation, do not take into account three-body corre¬ 
lations, and other methods must be developed. As dis¬ 
cussed above, strongly interacting Bose gases can reach 
an equilibrium state in the nondegenerate limit where the 
thermal wavelength is smaller compared to the interpar¬ 
ticle spacing. This regime should be accurately described 
by a virial expansion, and in this paper, we develop the 
virial expansion for the strongly interacting Bose gas in¬ 
cluding three-body corrections. 

This paper is structured as follows: in Sec. [TIJ we 
discuss the virial expansion and introduce the field- 
theoretical model for cold bosons that we shall use 
throughout the paper. Following this, we outline how 
three-body effects of the Bose gas are taken into ac¬ 
count in the virial expansion in a diagrammatic frame¬ 
work. Sectio n |III| presents the results of our work: we 
begin in Sec. |III A| by discussing our results for the first 
three virial coefficients. The second virial coefficient 
depends only on the scattering length and interpolates 
between a scattering-dominated limit for negative scat¬ 
tering length and a dimer-dominated limit for positive 
scattering length. The third virial coefficient depends 
not only on the scattering length but also sensitively on 
the three-body parameter k*. In particular, we find a 
rapid variation in the virial coefficient as the scattering 
length is tuned across the three-body thresholds where 
the system can support a single trimer bound state. In 
Sec. |IIIB[ we provide results for the two-body and three- 
body contact parameters, which describe various thermo¬ 
dynamic properties of the Bose gas and set the magni¬ 
tude of the asymptotic form of many response functions. 
Beside these thermodynamic quantities, we also develop 
the virial expansion for the full momentum distribution 
including three-body correlations, and the calculation of 
the momentum distribution, which constitutes the main 
result of our work, is presented in Sec. |III C| We discuss 
in detail the behavior of the momentum distribution as a 
function of both the scattering length and the three-body 
parameter. Our results are in excellent agreement with 
universal relations that govern the high-momentum tail 
of the momentum distribution, where we determine the 
contact parameters from the independent calculation pre¬ 
sented in Sec. |IIIB| Especially, we find very good agree¬ 
ment with a subleading log-periodic high-momentum tail 
which is a direct manifestation of Efimov physics. We 
also briefly make a comparison with a calculation of the 
Fermi momentum distribution, for which the virial ex¬ 
pansion provides only a quantitative correction but does 
not show new qualitative three-body effects. Our re¬ 
sults for the bosonic momentum distribution can be di¬ 
rectly compared with the recent experiment by Makotyn 
et al. ED, which is in the nondegenerate regime follow¬ 
ing a quench to the unitary limit. We provide a fit to the 
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FIG. 1. (a) Integral equation for the three-body scattering 
matrix, denoted by a gray box. Simple lines denote atom 
propagators, and double lines are dimer propagators, (b)-(m) 
Diagrams that contribute to the density and the momentum 
distribution up to third order in the fugacity. 


experimental data which contains only the temperature 
as a single fitting parameter. Our calculation provides 
an accurate quantitative description of the experimental 
results. Surprisingly, the best fit is found when the low¬ 
est Efimov trimer state is excluded from the calculation, 
indicating that the lowest trimer state is not occupied in 
the experiment ED- The paper is concluded with a sum¬ 
mary and conclusion in Sec.[V] There are two appendices 
which contain details of the calculation: Appendix [A] de¬ 
scribes the solution of the three-body scattering ampli¬ 
tude, and Appendix [B] contains intermediate results for 
diagrams that contribute to the virial expansion of the 
Bose gas. 


II. METHODS 


In this paper, we characterize the strongly interact¬ 
ing Bose gas in the normal phase by performing a virial 
expansion, allowing us to link few-body physics to the 
properties of an interacting many-body system in a sys¬ 
tematic way. In particular, in addition to the thermody¬ 
namic virial coefficients, we compute the full momentum 
distribution and establish Efimov signatures in this cor¬ 
relation function. The virial expansion is valid if the 
thermal wavelength is much smaller than the interparti¬ 
cle spacing (nX^ -C 1) and is therefore ideally suited to 
describe the current experimental work . Any op¬ 

erator can be expressed as a sum of separate traces over 
connected iV-particle sectors: 


Tr Oe-^n-^ 
Tre -/3(ff^A0 


£ z N tr N Oe-? H . 
N—0 


(i) 


In the nondegenerate regime nXj~ <C 1, the fugacity 
^ is a small parameter and the expansion in Eq. Q 
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can be truncated after the first few terms. By perform¬ 
ing the expansion up to third order in z, we fully in¬ 
clude three-body correlations. Previous early work on 
the virial coefficients of Bose gases is mum and more 
recent applications are [3TH13] . A previous numerical 
study of the virial coefficients was carried out by Be- 
daque and Rupak m, but they have not included all 
the relevant diagrams (see Fig. [lj. Note that the virial 
coefficients at unitarity can even be determined analyti¬ 
cally m ■ The virial expansion has also been successfully 
applied to Fermi gases; see Ref. |U] for a review. In the 
following, we provide numerical results for the first three 
virial coefficients (both at unitarity and at finite scatter¬ 
ing length) as well as the two-body and the three- body 
contacts. An essential new point of our analysis is that 
we develop the virial expansion for the full Green’s func¬ 
tion, which allows us to compute the momentum distri¬ 
bution. The momentum distribution exhibits a universal 
high-momentum tail that shows direct three-body Efi¬ 
mov correlations. We compare our calculations with the 
recent experiment m- 

The Lagrangian density of an interacting Bose gas with 
large scattering length takes the form mug 

C = $ (idt + <j> + yj -d^d - ^ 

~ ^^dU(f), (2) 


where cjr creates a boson and d, is an auxiliary dimer 
field. The bare coupling constants <?2 and <73 can be 
written as ^ ^ - rgw and g 3 = with 

m) ~ *.-e a »a cutoff r eg - 

ulator and A* is a renormalized three-body parame¬ 
ter For positive scattering length, there is a 

dimer bound state with energy Ejj — In addition, 

theory ([ 2 ]) has an infinite number of arbitrarily deep Efi¬ 
mov bound states. We avoid the Thomas collapse m 
of the zero range model by specifying a lowest Efimov 
trimer state, the energy of which we denote by Et■ The 
parameter A* is related to the wavenumber of the low¬ 
est Efimov trimer at unitarity via k* = 0.38A*. Each 
trimer branch hits the continuum of scattering states at 
the three-body threshold scattering length a_ and ter¬ 
minates at positive scattering length at the atom-dimer 
threshold a*. Our calculation is based on a diagram¬ 
matic representation developed for the virial coefficients 
of a Fermi gas [35], which we generalize to the Bose case. 
Crucially, we show that this method can be applied to any 
correlator and we identify the three-body effects on the 
momentum distribution. The starting point is an expan¬ 
sion in imaginary time, where it turns out that the full 
dependence on the fugacity z is encoded in the free prop¬ 
agator, which can be expanded in powers of z. In Fig. [T] 
we denote the jth-order contribution to the propagator 
by a line that is slashed j times. The zero-order con¬ 
tribution can only propagate forward in imaginary time, 
and backward-propagating lines contribute higher pow¬ 
ers of z. This provides a diagrammatic representation 
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FIG. 2. (Color online) (a) Second virial coefficient as a 
function of scattering length, (b) Third virial coefficient at 
unitarity as a function of the three-body parameter. The 
dashed red and dotted green lines in (a) and (b) denote the 
analytical asymptotic results [411 132] . (c) Third virial coef¬ 
ficient as a function of scattering length for fixed (along the 
arrow) k*A t = 0.5,1,2,3, 5. (d) Bound-state spectrum. Blue 
solid lines denote trimer branches, and the red dashed line is 
the dimer bound-state energy, (e) Third virial coefficient at 
small and positive scattering length a < At. Dashed lines 
denote the asymptotic result based on the atom-dimer Beth- 
Uhlenbeck formula 


for the virial expansion of any correlation function. The 
Feynman diagrams constructed in this way contain sub¬ 
diagrams that involve repeated scattering of either two 
or three particles. These are the scattering matrices of 
the few-body problem, and we denote them by a double 
line and a box, respectively. 


III. RESULTS 


A. Virial coefficients 


We begin by computing the virial expansion of the den¬ 
sity, which is related to the virial coefficients via [35] 


X^n = b\Z + 2& 2 ~ 2 + 3& 3 z 3 + 0(z 4 ). 


(3) 


The diagrams contributing to the d ensity up to 0(z 3 ) 
are shown in Fig. lj with Fig. |l(b)| contributing to 61 , 
Figs. |l(cj| and 1(e) contributing to 62 , and all other dia¬ 
grams contributing to the third virial coefficient 63 . Fig- 
l(b)|l(d) give the standard result b^ = j ~ 5 / 2 
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FIG. 3. (Color online) (a) Two-body contact as a function of 
temperature. Three-body contact at unitarity as a function of 
(b) temperature and (c) the three-body parameter. The lines 
in (a) and (b) correspond to (from bottom to top) K*/fc n = 
1,2,3,4, and 5. The green lines in (a) and (b) denote the 
asymptotic high-temperature results. The red and green lines 
in (c) denote the asymptotic results as stated in the main text. 


FIG. 4. (Color online) Third-order contribution n 3 {Q to the 
momentum distribution at unitarity with ft* At = 3 as a func¬ 
tion of dimensionless momentum £ = q At- Points denote 
numerical results; the solid blue line is a guide to the eye. 
The green and dot-dashed orange lines denote the asymp¬ 
totic Tan relations (Ffl obtained with e~ l 3 ET C 2,3 = 825.2 and 
e~ pBT c 3 =48.1. 


for the noninteracting Bose gas, while Fig. |l(e)| yields 
the Beth-Uhlenbeck interaction correction m 09]. Fig¬ 
ures |l(i)fl(ni) have to be evaluated numerically, and 
the three-body amplitude in the integrand is obtained 
by solving the Skorniakov-Ter-Martirosian integral equa¬ 
tion ESI shown in Fig. m Extensive details of our cal¬ 
culation are relegated to the appendices. Figure [2] shows 
our results for the second and third virial coefficients. 
Both the second [Fig. [2j[a)] and third virial coefficients 
[Fig. |2jb)] agree with the analytical results HD EH- pro- 
viding an independent check of our calculations. Fig¬ 
ure E c ) shows the third virial coefficient as a function of 
scattering length for various values of ft* A t- We rescale 
the scattering length by a_ for each ft* and plot the re¬ 
duced virial coefficient e~ t3E b 3 , where E is the lowest 
bound-state energy (E = E T for a -1 < aj 1 , E = E D for 
a -1 > a” 1 , and E = 0 if there is no bound state). In 
the scattering-dominated limit ft*A t 1 , the virial co¬ 
efficient increases smoothly as the scattering length is 
tuned across the unitary limit and then decreases for 
a -1 > a~ D assuming negative values. In the trimer 
limit k*At 1 , the reduced virial coefficient is very 
small if no Efimov bound state is contained in the spec¬ 
trum and rapidly jumps to the trimer-dominated result 
e -/3 Er b 3 — 3\/3 [39] 05] otherwise. Figure [TJe) magnifies 
63 at small scattering length a -1 > a~ 1 for which only 
the dimer bound state exists. In this limit, we find that 
the virial coefficient is well described by an effective two- 
body Beth-Uhlenbeck formula for the scattering of atoms 


and dimers 09] EE] 



where a a a («,«*) is the atom-dimer scattering length [J|. 
Note that the zero-range description ([ 2 ]) holds if the 
scattering length and the thermal wavelength are large 
compared to the effective range: a, A t > i v dw (how¬ 
ever, a and A t can have arbitrary ratios). Since the 
three-body parameter is found to be universal with 
ft* « §.2/l vdw [52][S3], we expect that in the scattering- 
dominated limit, where UbT is larger than the energy of 
the lowest Efimov trimer, ft*A t 1, there are effective 
range corrections to the results of this paper. 


B. Two- and three-body contacts 


There exists a large set of universal relations that de¬ 
scribe thermodynamic quantities and the short-time and 
short-distance structure of a quantum gas [53H59] . The 
magnitude of these relations is set by the two-body con¬ 
tact C 2 and the three-body contact C 3 , which are re¬ 
lated to the energy of the system by the adiabatic theo¬ 
rems [53H59] : 


C2 = —87TTO 


dE 

da - 1 


mn * dE 
2 5ft* _i 


T j >2 

(5) 

^ V" i 

At j>3 

(6) 
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FIG. 5. (Color online) Full momentum distribution n(() up to third order in the fugacity with k»At = 3 for three different 
scattering lengths A t/cl = —3.16,0, and 3.16. Fugacities have been chosen such that z 3 e l3E = 0.6. The notation is the same as 
in Fig. [4] 





FIG. 6. (Color online) Full momentum distribution n(() up to third order in the fugacity at unitarity with k»A t = 1, 3, and 5. 
Fugacities have been chosen such that z 3 e ,3E = 0.6. The notation is the same as in Fig. Ul 


where c 2 j = 16n 2 dbj /8(Xt/ a) and c 3 j = 7 TK*dbj /8k*, 
and we abbreviate 03,3 = C 3 . In Fig. [3] we plot the 
intensive contact parameters C 2 /Nk n [Fig. ia)] and 
C 3 /Nk 2 n [Fig. §b)} as a function of rescaled tempera¬ 
ture T/T n , where k n = (6ir 2 n) 1//3 and T n = k 2 /2m. For 
T T„, the contacts approach the asymptotic results 
C 2 /Nk n = 64T„/3T and C 3 /Nk 2 = is 0 T 2 /VSn 2 T 2 
(green lines without dots in Fig. |) As the temperature 
is lowered, both C 2 and C 3 strongly increase and then 
approach a constant value at low temperature. In ad¬ 
dition, we show the three-body contact as a function of 
k*A t in Fig- [3|c). The numerical results agree with the 
analytical limits C 3 = 3v / 3(K*AT) 2 e /3BT for large k*A t 
and C 3 = 3\/3so for small k*A t Ea¬ 


ch Momentum distribution 

A particular universal relation connects the two- and 
three-body contacts to the high-momentum tail of the 
momentum distribution [S3EZU5E!: 

n( q ) = ^ + ^ F ( q ) + °(l/q e ), (7) 

where C 2/3 = C 2 / 3 /V is the intensive contact density and 
F(q) is a log-periodic function of the momentum given 
by -F(g) = Asin( 2 soln( 7 /ft* + 20), where A = 89.26260 
and 0 = —0.669064. Indeed, a comparison of the exper¬ 
iment m with the asymptotic result ([7]) using a scal¬ 
ing ansatz for the contact parameters suggests that the 
observed momentum distribution is consistent with the 
three-body tail m■ Note that the momentum distribu¬ 
tion can be computed for a single trimer [5BJED]- Here, 
we provide a full calculation of the momentum distribu¬ 
tion at finite density to third order in the fugacity. 
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The momentum distribution n(q) is defined as the zero¬ 
time limit G(0~,q) of the imaginary-time propagator. 
Formally, the diagrams that contribute to n(q) up to 
third order in z are the same as for the density in Fig. [T] 
although the black dot now denotes a momentum inser¬ 
tion q into the diagram, changing the structure of the 
calculation fundamentally. Figure [4] shows our result for 
the 0(z 3 ) part n 3 (q) of the momentum distribution at 
unitarity with k*At = 3, where we parametrize 

n(q) = e^ ET [ni(q)z + n 2 {q)z 2 + n 3 {q)z 3 + 0(z 4 )]. (8) 

Figure |4ja) shows n 3 (q) as a function of dimension¬ 
less momentum = qXr, Fig. |4jb) amplifies the high- 
momentum tail, and Fig. [djc) shows the momentum dis¬ 
tribution with a subtracted leading-order two-body term. 
Our results are in excellent agreement with the leading- 
order (red dashed line) and next-to-leading order (green 
line) Tan relations where the 0{z 3 ) contact param¬ 
eters are extracted from the virial coefficients in an in¬ 
dependent calculation, providing a strong check of our 
results. We remark that the detailed comparison of the 
numerical high-momentum tail with the independently 
obtained contact parameters over more than three orders 
of magnitude in momentum poses exceptional constraints 
on the numerical implementation, requiring a relative er¬ 
ror of 10 -9 in our simulations. Note that the onset of the 
asymptotic tail is set by k K = \f2mEr rather than k n , 
which marks the asymptotic regime in deeply degenerate 
Fermi gases m- 

To illustrate the behavior of the momentum distribu¬ 
tion as a function of scattering length and the three-body 
parameter, we show the full momentum distribution in¬ 
cluding all terms up to third order in the fugacity in 
Fig -0 for fixed «*A t = 3 and three different scatter¬ 
ing lengths cl/\t — —3.16,0, and 3.16. In Fig. [ 6 ] the 
momentum distribution at unitarity is shown for three 
values of the three-body parameter, k*Xt = 1,3, and 5. 
In both Figs. [5] and [ 6 ] the fugacity was chosen such that 
z 3 e. 3E = 0.6, where E > 0 again denotes the deepest 
bound-state energy. All momentum distributions agree 
with the contact parameters at unitarity or finite scat¬ 
tering length as obtained from the adiabatic theorems 
using the virial coefficients, which were computed in an 
independent calculation. Figures [5 (a) | and |5(c)| illustrate 
that, compared to the unitary case, the presence of an¬ 
other scale ( 1 /a ^ 0 ) delays the saturation to the univer¬ 
sal high-momentum behavior ([T]) by almost an order of 
magnitude in momentum. 

To conclude this section, we briefly make a comparison 
with the fermionic case, which can be computed along the 
lines outlined in the appendices. For an approach that 
calculates the momentum distribution via the spectral 
function, see [B2]. Figure [T] shows the momentum dis¬ 
tribution n|(C) of spin-up particles in a two-component 
Fermi gas with equal population in both spin species, cal¬ 
culated to third order in the fugacity. The three different 
scattering lengths are \t/ a = — 1 , 0 , 1 , and the fugacities 
were chosen such that z 2 e^ E = 0.15. While at first glance 



C — qX t 

FIG. 7. (Color online) Momentum distribution of spin-up 
particles in a population-balanced two-component Fermi gas 
as a function of dimensionless momentum £ = qXr including 
terms up to third order in the fugacity z. The scattering 
lengths are A r/a = —1 (green), Xt/cl = 0 (red) and A r/a, = 1 
(blue). Fugacities were chosen such that z 2 e 3E = 0.15, where 
E = Ed for a > 0 and E = 0 otherwise. Dot-dashed lines 
denote the corresponding values of the contact as obtained 
from an independent calculation for the virial coefficients. 



FIG. 8. (Color online) Momentum distribution of a unitary 
Bose gas. k = q/k n denotes the dimensionless momentum. 
The solid red and blue lines are the result of the trap-averaged 
virial expansion with fugacities z\ = 0.5 and z 2 = 0.4, re¬ 
spectively, corresponding to the experimental data m at 
(m) = 5.5-10 12 cm -3 (orange line) and ( 712 ) = 1.6-10 12 cm -3 
(green line). The dashed red and blue lines are the corre¬ 
sponding 0(z 2 ) results. All momentum distributions are nor¬ 
malized to unity J d 3 Kn(n)/(2n) 3 = 1. 


the momentum distributions look qualitatively the same 
as the bosonic ones, the fermionic system shows no os- 
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cillatory behavior in its high-momentum tail, as is to be 
expected since the three-body contact is without analog 
in fermionic gases with equal masses. 


IV. COMPARISON WITH EXPERIMENTS 


We compare our results with the experiment by Mako- 
tyn et al. EB, which measures the momentum distri¬ 
bution following a quench to the unitary limit for two 
different initial densities (ni) = 5.5 • 10 12 cm -3 and 
(n 2 ) = 1.6 • 10 12 cm -3 . Following Ref. (43], we assume 
a constant phase space density nX^ and average our re¬ 
sults over a Thomas-Fermi density profile. Since our cal¬ 
culations are performed at fixed ft*A t, we keep ft*/fc n 
fixed at its value at the trap center, which neglects loga¬ 
rithmic corrections to the fugacity. The experiment was 
performed using 85 Rb, which has a three-body parame¬ 
ter ft* = 38(1) pmT 1 (631 . Remarkably, our results are in 
good agreement with the experiment if and only if we ex¬ 
clude the lowest trimer branch from our calculation and 
set ft* = ft*/22.7, indicating that on the time scales of the 
experiment |31| . the lowest Efimov trimer branch is not 
occupied. In this case, the virial expansion agrees well 
with the experimental data with Z\ = 0.5 and z 2 = 0.4. 
The small values of the fugacity justify the use of the 
virial expansion. The results are shown in Fig. [ 8 ] For 
comparison, we include a fit up to second order in z as a 
dashed line. 


V. RESULTS AND CONCLUSIONS 


In summary, we have characterized the strongly inter¬ 
acting Bose gas in the normal phase by computing the 
first three virial coefficients, the two-body and the three- 
body contacts, and the momentum distribution. The re¬ 
sults for the momentum distribution are in good quan¬ 
titative agreement with the recent experiment |31i . Sur¬ 
prisingly, the fit indicates that the lowest Efimov trimer 
state is not populated in the experiment. Our work opens 
the possibility for the spectroscopy of Efimov states at 
unitarity. 
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Appendix A: STM equation 


The 0(z 3 ) of the virial expansion relates the density 
and the momentum distribution to the vacuum three- 
body T 3 matrix, and we summarize some basic proper¬ 
ties in this section. The T 3 matrix is obtained by solv¬ 
ing the Skorniakov-Ter-Martirosian (STM) integral equa¬ 
tion [50j . which is shown diagrammatically in Fig. 1(a). 
The three-body matrix depends on the total energy s 
and the total momentum P of the incident atom-dimer 
state as well as the energy and momentum ( E in , k in ) and 
(E out ,k out ) of the ingoing and outgoing atom line. It 
turns out that our calculations require only the on-shell 
T 3 matrix where the external atom energy is equal to the 
kinetic energy of a free atom, i.e., E in = = k 2 /2m 

and E ou t = e p = p 2 /2m. Using the Feynman rules of the 
Lagrangian (2) in the main text, the integral equation for 
the on-shell T 3 matrix with total energy s and incoming 
and outgoing atom momenta p and k reads 


Tz(s, £ P , £k, P, p, k) 


1 

S £p £k £P—p—k 




x T 2 (s-£q- T 3 (s,£ q ,£ k ,P,q, k), 

(Al) 


where we have defined f k = f d 3 k/(2 tt) 3 as a shorthand 
for the momentum integrals. T 2 is the two-particle scat¬ 
tering matrix in vacuum: 


T 2 (s) 


8-7T 1 

mA — y/—ms 


(A2) 


Due to the Galilean invariance of the vacuum state, we 
have the following important transformation property: 


l3(s,£ p ,£ k ,P,p,k) = h 





(A3) 


where £ 3 ( 5 , q,q') = T 3 (s, £ q , £ q ', 0, q, q') denotes the T 3 
matrix in the center-of-mass frame. The STM equation 
can be decoupled in angular momentum channels using 
the angular decomposition 

t 3 (s,p,k) = 5^(2 1 + l)Pi(cos6)t%\s,p,k) (A4) 

l 

1 

t%\s,p,k) = -J dcos 9 Pz(cos 0)i 3 (s, p, k), (A5) 


where Pi (cos 9) are the Legendre polynomials and cos 9 = 
p • k, where the hat denotes a unit vector. This yields a 
set of decoupled STM equations for the different angular 











momentum channels: 


4 Z) ( S >-P> fc ) = 


pk^ 1 \pk 


k 2 


S _P_ _ 

m m 


H(A) ■ 


- / dq- 


- — \ —ms - 


!<? 2 


~-t^\ s ,q,k) 


m 


—Qi - 


pq 


m 


PQ 


s _^_^ 
m m 


H (A) ' 
- m-p-Sm 


■ (A 6 ) 


Here, Qi(z) denotes a Legendre function of the second 
kind [5J. As discussed in the main text (see also [!]), 
we renormalize the theory in order to eliminate the bare 
coupling constants in favor of the log-periodic cutoff func¬ 
tion H{ A). The numerical solution of the integral equa¬ 
tion ( |A 6 | ) in the l = 0 channel yields an Efimov spectrum 
with infinitely many three-body bound states and a dis¬ 
crete scaling symmetry as described in the Introduction 
of the main text. 


distribution n(q), the external lines have a momentum q 
which is not integrated over. 

Like for fermions @Sj, the j th order in the virial ex¬ 
pansion of a correlator can be constructed systematically 
by constructing the Feynman diagrams with a maximum 
number of j backward-propagating lines. These diagrams 
involve an infinite number of repeated scattering of two or 
three forward-propagating lines, which, when resummed, 
yield the vacuum two-body T 2 and three-body X 3 matri¬ 
ces. The virial expansion of the density gives the virial 
coefficients and, using the adiabatic theorems, the two- 
and three-body contact parameters as discussed in the 
main text. Our calculation, which is set up for the full 
Green’s function, allows us to go beyond these results to 
compute the full momentum distribution, which contains 
universal Efimov correlations. In the following, we dis¬ 
cuss the derivation of the diagrammatic contributions to 
the momentum distribution up to third order in z and 
state the results for all diagrams. 


Appendix B: The momentum distribution 


2. One-body diagrams 


1. Expansion of the Green’s function 

In the following, we will illustrate the calculation for 
the momentum distribution n(q). The momentum distri¬ 
bution is defined as the zero-time limit of the imaginary- 
time Green’s function: 

n(q) = — lim G(r, q). (Bl) 

T —^0 

Tire total density is defined as 

n = [ n( q) = - lim_ f G(r, q), (B2) 

J q T_>0 J q 

which implies that the results of the diagrams for the 
density can be obtained via integration of the results for 
the momentum distribution over all momenta. The non¬ 
interacting Green’s function Go can be expanded in the 
fugacity by expanding the Bose-Einstein distribution in 
powers of z: 


G 0 (r,q) = -e- T ^-^ 


0(r) + n B (Eq - p) 


= q). 


(B3) 


3 =0 


Diagrammatically, the 0(z n ) contribution g\P (t. q) in 
the expansion of the free Green’s function is denoted by 
a line that is slashed n times. Since the zero-order contri¬ 
bution is purely retarded, all backward-running lines con¬ 
tribute at least one power of z to a diagram. The terms 
that define the right-hand side of Eq. (B31 up to third 


order in the fugacity yield the one-body diagrams shown 
in Figs. l(b)-l(d) of the main text. For the momentum 


The contribution to a one-body diagram that is slashed 
^ times can be directly inferred from Eq. (B31: 


lim G^(T,q) = -e-^, (B4) 

T —^0 


where Fig. 1(b) has i = 1, Fig. 1(c) corresponds to 
1= 2, and Fig. 1(d) corresponds to £ = 3. 


3. Two-body diagrams 

We will consider the 0(z 2 ) two-body diagram in some 
detail, and quote the results for the remaining diagrams 
for completeness. We denote the imaginary-time differ¬ 
ence between time zero and the scattering vertex by t±, 
while the scattering vertex shall act over a time difference 
t 2 ■ With these conventions, Fig. 1(e) reads: 


= —z 2 dt\dt 2 Go°\ti, qjGg 1 ) (—t 2 , k) 

x e 2/it2 T 2 (f 2 ,q + k)Gp 1) (-fi - 1 2 ,q) 

= z 2 f dt x dt 2 f T 2 (t 2 , q + k)e- tl ( e - +Ek > 

J A J q 

x (B5) 

where A = {(ti,f 2 ) : 0 < t x + t 2 < /3,0 < ti,0 < t\}. 
Using the generalized convolution theorem for Laplace 
transforms m ns], we can write this as the inverse 
Laplace transform of the Laplace transforms of the three 
functions in the integrand: 



-/ 3 S /' T 2 (s,k + q) 

J k (s-e k -e q ) 2 ' 


(B 6 ) 
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The integration f s = J BW ds/2ni denotes a Bromwich 
contour in the complex energy plane, which runs (parallel 
to the imaginary axis) to the left of the dimer pole and 
the branch cut of the T 2 matrix. We now substitute 
p = q + k and s — > s — £ p , which allows us to eliminate 
the angle integrations from the momentum integration 
and to transform the above expression into 


metry factor of 1/2, while Figs. l(j)-l(m) share the 
same factor of 1. In addition, the lack of the integra¬ 
tion over the external momentum considerably compli¬ 
cates the evaluation of the diagrams both numerically 
and analytically compared to the density. Eliminating 
the convolution integral over the imaginary times in fa¬ 
vor of a single complex integration, we find for Fig. l(i) 



(2tt)2 


= 2 2 / e-^T 2 (s) 

S 

d P , , 2 


-l £ 
e 2 fc p 


s Jo 


jp ( s d 2 e q-p e q)“ 

2 e-P s T 2 {s)p 2 e-P^ 

( 2 - s] 2 - ^ + s] 

y 4m J ^ ^m °) 2 m \m ' a ) 

(B7) 


The remaining two-body diagrams can be evaluated in a 
similar way, and we list only the results: 




= ze 



(B 8 ) 



= 2 *11 dp 


p2 e P(s+ + 6 a m ) T 2 (s) 


s JO 


fs-Pi') _ ( s + P 2 -') + 

4m J 9m ( T 4raj ' 


i6 r 

81 m 2 


(B9) 


= zr e 


I k.p 


T„ 2 (s-£ q ,k + p) 


Pi (S; £q; ^q, k 


p + q,q,q) 


(s £k £ p £ q )“ 


(BIO) 


where k and p denote the loop momenta of the two 
backward-running atom propagators that are connected 
to the dimer fields. In the case of the three-particle di¬ 
agrams, the Bromwich contour runs to the left of the 
trimer pole which we have chosen to be the physically 
deepest one. This amounts to an implicit subtraction of 
all deeper trimer poles. We now substitute 


s = s — 


£q+p+k 


-±1, 


— sl.i sla 


-P3 


p - q 

k-q 


- 3 


2 0 _ llr _ Ip, 

3 q 3 K 3 P 

Ik - In 

2 K 2C 


(Bll) 


(B12) 


4. Three-body diagrams 

The three-body diagrams differ in various ways. First 
and foremost, note that Fig. l(i) gets an overall sym- 


where the determinant of the transformation contributes 
a factor 3 3 and I 3 denotes the identity matrix in three 
dimensions. Using Eq. (A31 and relabeling s' —> s,k' —► 
k,p' —> p, the contribution of Fig. l(i) can be expressed 
in terms of the center-of-mass amplitudes as 



_3V f Pa f „-3a En _ r 2 2 (g-|gp)^3(^P-P) 


' k,p 


(s - j£p - 2£ k ) 2 


3 2 z 3 m 3 
4g/3(27r) 3 


X ](2 ^ + 1 ) 


dpp- 


T 2 I a _ Ip1 

±0 1 S 4m 


) t$\s,p,p) 


e -^U 2 +P 2 ) s inh 




( 3 pqP\ 


-ms + 


\ m 


■ 


(B 13 ) 


For the last line, we have solved the elementary integral over k and used the fact that 43 (s, p, p) = -Pz(P'P)4^( s jIM') 

contains only forward scattering, such that all the Legendre polynomials equal to 1. We have then integrated over 
the angles of p. We will now list the resulting expressions for the other diagrams, which are obtained in an analogous 
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fashion: 


3 1 2 z 3 m 


(27t) 4 /3 t—j 


p pOO pOO 

(:21 + 1 ) / e~ l3s / dp dkpkt^\s,k,k) 
Js Jo Jo 


e~i^( k ~p'> 2 — e ~§^(p+ k ) 2 

(«- ^ 2 ) 


3k 2 , (P~2q) 2 

4m ' 4m 


3k 2 , (p+2q) 2 

4m ' 4m 



(B14) 



2m 3 3 2 2 3 
fdq (2tt) 4 


X ](2 ^ + 1 ) 


— / dfcT2 ( s ~ 4^ fc2 ) T2 ( S_ 


/o P Jo 


e -|^(p- 9) 2 _ e -|^(P +?) 2 




4m 

p 2 k 2 
m m 


4°( S >P, fc), 


where Qz(^) = —dQi(z)/dz is the solution of the angle integration. Furthermore, we have: 



2m 3 3 2 z 3 


fiq (2tt) 4 ^ 


5^(2 1 + 1) 


r ± r 

1 0 P Jo 


3 

4m 


n -t>e I -r l dkT ^ s _ OPtf )T 2 ( a _ ^ 

p 2 k 2 

s - 

m m 


e -J i :( p - 9) 2 _ e -|^( P +?) 2 


Qi 


pk 


4 Z) ( s ’ fc >p)- 


For the last diagram, the integration over the angles can only be partially performed analytically: 



3V 2m f J^e-P s r dkk 2 r dpP 2 


(27r) 4 3 pq J 2iri 


dcosd 


(B15) 


(B16) 


e-^ {k +p +9 ) T 2 (s-—k 2 )T 2 (s-—p 2 )Y J ^l + m\s,k,p) 


Am 


Pi(cos9)e 


-M kpcos6 sinh \A 2 + P 2 + 2 pk cos 9 


(s — — — — - ^COS0) 2 

v m m m > 


Vk 2 + p 2 + 2 pk cos 9 


(BIT) 


5. Numerical evaluation 


During the numerical evaluation, we have retained an¬ 
gular momenta up to l = 10, finding no change when ically over the loop momenta and the complex energies 

including higher harmonics into the calculation. In each using the numerical solution of the STM equations, 

angular momentum channel, the STM equation (A6) 


r 


is solved at different complex energies that lie on the 
Bromwich contour. Afterwards, all Feynman diagrams 
discussed in the previous sections are integrated numer- 


[1] M. Suvakov and V. Dmitrasinovic, Phys. Rev. Lett. 110, 
114301 (2013). 

[2] V. Efimov, Physics Letters B 33, 563 (1970). 

[3] V. Efimov, Sov. J. Nucl. Phys. 12, 589 (1971). 

[41 E. Braaten and H.-W. Hammer, Physics Reports 428, 
259 (2006). 

[5] T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, 
C. Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, 


H. C. Nagerl, and R. Grimm, Nature 440, 315 (2006). 

[6] S. Knoop, F. Ferlaino, M. Mark, M. Berninger, 
H. Schobel, H. C. Nagerl, and R. Grimm, Nat. Phys. 
5, 227 (2009). 

[7] M. Zaccanti, B. Deissler, C. D’Errico, M. Fattori, 
M. Jona-Lasinio, S. Muller, G. Roati, M. Inguscio, and 
G. Modugno, Nat. Phys. 5, 586 (2009). 






































11 


[8] N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, 
Phys. Rev. Lett. 103, 163202 (2009). 

[9] S. E. Pollack, D. Dries, and R. G. Hulet, Science 326, 
1683 (2009). 

[10] T. B. Ottenstein, T. Lompe, M. Kohnen, A. N. Wenz, 
and S. Jochim, Pliys. Rev. Lett. 101, 203202 (2008). 

[11] J. H. Huckans, J. R. Williams, E. L. Hazlett, R. W. 
Stites, and K. M. O’Hara, Phys. Rev. Lett. 102, 165302 
(2009). 

[12] J. R. Williams, E. L. Hazlett, J. H. Huckans, R. W. 
Stites, Y. Zhang, and K. M. O’Hara, Phys. Rev. Lett. 
103, 130404 (2009). 

[13] T. Lompe, T. B. Ottenstein, F. Serwane, A. N. Wenz, 
G. Ziirn, and S. Jochim, Science 330, 940 (2010). 

[14] S. Nakajima, M. Horikoshi, T. Mukaiyama, P. Naidon, 
and M. Ueda, Phys. Rev. Lett. 106, 143201 (2011). 

[15] G. Barontini, C. Weber, F. Rabatti, J. Catani, G. Thal- 
hammer, M. Inguscio, and F. Minardi, Phys. Rev. Lett. 
103, 043201 (2009). 

[16] S.-K. Tung, K. Jimenez-Garcia, J. Johansen, C. V. 
Parker, and C. Chin, Phys. Rev. Lett. 113, 240402 
(2014). 

[17] R. Pires, J. Ulmanis, S. Hafner, M. Repp, A. Arias, E. D. 
Kuhnle, and M. Weidemiiller, Phys. Rev. Lett. 112, 
250404 (2014). 

[18] M. Kunitski, S. Zeller, J. Voigtsberger, A. Kalinin, 
L. P. H. Schmidt, M. Schoffler, A. Czasch, W. Schollkopf, 
R. E. Grisenti, T. Jahnke, D. Blume, and R. Dorner, Sci¬ 
ence 348, 551 (2015). 

[19] Y. Nishida, S. Moroz, and D. T. Son, Phys. Rev. Lett. 
110, 235301 (2013). 

[20] Y. Nishida, Y. Kato, and C. D. Batista, Nat. Phys. 9, 
93 (2013). 

[21] V. Efimov, Nuclear Physics A 210, 157 (1973). 

[22] D. S. Petrov, Phys. Rev. A 67, 010703 (2003). 

[23] N. Navon, S. Piatecki, K. Gunter, B. Rem, T. C. Nguyen, 
F. Chevy, W. Krauth, and C. Salomon, Phys. Rev. Lett. 
107, 135301 (2011). 

[24] L.-C. Ha, C.-L. Hung, X. Zhang, U. Eismann, S.-K. 
Tung, and C. Chin, Phys. Rev. Lett. 110, 145302 (2013). 

[25] P. O. Fedichev, M. W. Reynolds, and G. V. Shlyapnikov, 
Phys. Rev. Lett. 77, 2921 (1996). 

[26] E. Nielsen and J. H. Macek, Phys. Rev. Lett. 83, 1566 
(1999). 

[27] B. D. Esry, C. H. Greene, and J. P. Burke, Phys. Rev. 
Lett. 83, 1751 (1999). 

[28] J. P. D’lncao, H. Suno, and B. D. Esry, Phys. Rev. Lett. 
93, 123201 (2004). 

[29] R. J. Fletcher, A. L. Gaunt, N. Navon, R. P. Smith, and 
Z. Hadzibabic, Phys. Rev. Lett. Ill, 125303 (2013). 

[30] B. S. Rem, A. T. Grier, I. Ferrier-Barbut, U. Eismann, 
T. Langen, N. Navon, L. Khaykovich, F. Werner, D. S. 
Petrov, F. Chevy, and C. Salomon, Phys. Rev. Lett. 110, 
163202 (2013). 

[31] P. Makotyn, C. E. Klauss, D. L. Goldberger, E. A. Cor¬ 
nell, and D. S. Jin, Nat. Phys. 10, 116 (2014). 

[32] U. Eismann, L. Khaykovich, S. Laurent, I. Ferrier- 
Barbut, B. S. Rem, A. T. Grier, M. Dclehaye, F. Chevy, 
C. Salomon, L.-C. Ha, and C. Chin, arXiv:1505.04523 
(2015). 

[33] H. T. C. Stoof, Phys. Rev. A 49, 3824 (1994). 


[34] E. J. Mueller and G. Baym, Phys. Rev. A 62, 053605 

( 2000 ). 

[35] G. S. Jeon, L. Yin, S. W. Rhee, and D. J. Thouless, 
Phys. Rev. A 66, 011603 (2002). 

[36] W. Li and T.-L. Ho, Phys. Rev. Lett. 108, 195301 (2012). 

[37] S. Piatecki and W. Krauth, Nat. Commun. 5, 3503 
(2014). 

[38] W. van Dijk, C. Lobo, A. MacDonald, and R. K. 
Bhaduri, arXiv:1412:5112 (2014). 

[39] A. Pais and G. E. Uhlenbeck, Phys. Rev. 116, 250 (1959). 

[40] R. Dashen, S.-K. Ma, and H. J. Bernstein, Phys. Rev. 
187, 345 (1969). 

[41] P. F. Bedaque and G. Rupak, Phys. Rev. B 67, 174513 
(2003). 

[42] Y. Castin and F. Werner, Canadian Journal of Physics 
91, 382 (2013). 

[431 S. Laurent, X. Leyronas, and F. Chevy, Phys. Rev. Lett. 
113, 220601 (2014). 

[44] X.-J. Liu, Physics Reports 524, 37 (2013). 

[45] P. F. Bedaque, H.-W. Hammer, and U. van Kolck, Phys. 
Rev. Lett. 82, 463 (1999). 

[46] P. Bedaque, H.-W. Hammer, and U. van Kolck, Nuclear 
Physics A 646, 444 (1999). 

[47] L. H. Thomas, Phys. Rev. 47, 903 (1935). 

[48] X. Leyronas, Phys. Rev. A 84, 053633 (2011). 

[49] E. Beth and G. E. Uhlenbeck, Physica 4, 915 (1937). 

[50] G. Skorniakov and K. Ter-Martirosian, Sov. Phys. JETP 
4, 648 (1957). 

[51] An analogous result holds for the fermionic case with 
a a d = 1.18a, see V. Ngampruetikorn, M. M. Parish, and 
J. Levinsen, Phys. Rev. A 91, 013606 (2015). 

[52] R. Schmidt, S. P. Rath, and W. Zwerger, EPJ B 85, 1 

( 2012 ). 

[53] J. Wang, J. P. D’lncao, B. D. Esry, and C. H. Greene, 
Phys. Rev. Lett. 108, 263001 (2012). 

[54] S. Tan, Annals of Physics 323, 2952 (2008). 

[55] E. Braaten, D. Kang, and L. Platter, Phys. Rev. Lett. 
106, 153005 (2011). 

[56] F. Werner and Y. Castin, Phys. Rev. A 86, 053633 

( 2012 ). 

[57] S. Tan, Annals of Physics 323, 2971 (2008). 

[58] Y. Castin and F. Werner, Phys. Rev. A 83, 063614 

( 2011 ). 

[59] D. H. Smith, E. Braaten, D. Kang, and L. Platter, Phys. 
Rev. Lett. 112, 110402 (2014). 

[60] F. F. Bellotti, T. Frederico, M. T. Yamashita, D. V. Fe¬ 
dorov, A. S. Jensen, and N. T. Zinner, Phys. Rev. A 87, 
013610 (2013). 

[61] J. T. Stewart, J. P. Gaebler, T. E. Drake, and D. S. Jin, 
Phys. Rev. Lett. 104, 235301 (2010). 

[62] M. Sun and X. Leyronas, Phys. Rev. A 92, 053611 (2015) 

[63] R. J. Wild, P. Makotyn, J. M. Pino, E. A. Cornell, and 
D. S. Jin, Phys. Rev. Lett. 108, 145305 (2012). 

[64] W. Magnus, R. P. Soni, and F. Oberhettinger, Formulas 
and theorems for the special functions of mathematical 
physics (Springer, 1966). 

[65] G. Doetsch, Introduction to the Theory and Application 
of the Laplace Transformation (Springer Science & Busi¬ 
ness Media, 2012). 


